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ABSTRACT 


A consistent linearization is provided for the element-independent corotational formulation, 
providing the proper first and second variation of the strain energy. As a result, the warping 
problem that has plagued flat elements has been overcome, with beneficial effects carried 
over to linear solutions. True Newton quadratic convergence has been restored to the 
Structural Analysis of General Shells (STAGS) code for conservative loading using the 
full corotational implementation. Some implications for general finite element analysis 
are discussed, including what effect the automatic frame invariance provided by this work 
might have on the development of new, improved elements. 
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1 . Introduction 


In this report, the problem of element warping and slow, irregular convergence to nonlinear 
solution states will be addressed in the context of the element-independent corotational 
theory [7] now implemented in STAGS. We shall show that some of these results have 
implications beyond the original confines of the theory, and beyond the SH410 elements 
that were used in this investigation. Section 2 addresses the problems encountered during 
a linear analysis related to warping, and Section 3 addresses questions of convergence to 
nonlinear solutions. 


1 


2. The Problem of Warping 

Many cases have surfaced where the normally reliable workhorse SH410 (410 for short) 
shell element has failed to produce a converged solution even for very fine grids. Perhaps 
the worst case of all is the large class of problems involving a doubly-curved surface, or in 
fact any curved surface for which the nodes do not coincide with the surface. Even though 
the finite element codes calculate the actual nodal coordinates, only the projected nodal 
positions on the plane formed by the element can be used in the element calculations: the 
410 is a flat plate element. For the vast majority of cylinder, annulus, and even the sphere 
problem with the usual polar grid, the element plate sections properly tile the surface and 
difficulties are avoided. Unfortunately, many problems exist where it is impossible to tile 
the surface exactly using flat elements, especially when grid refinement is necessary near 
discontinuities, or when design features spoil the simple cylinder or sphere geometry. An 
example of this is a curved composite plate with a cutout, where the cutout plays the role 
of spoiler. If these plates have significant curvature, unacceptable inaccuracies could be 
expected. 

Approximately a year ago, when the four and nine node assumed natural strain (ANS) 
elements were being developed [5] , the often-used pinched hemisphere problem was chosen 
to check the elements out. Two grids were tried. The first was the usual polar grid with 
the single row of triangles near the apex, with the remainder as quadrilaterals exactly 
tiling the surface along latitude and longitude lines. The second was the grid shown in 
Figure 1, where a quarter of the sphere was modeled with symmetry along the sides and 
a free edge at the base. Note that this grid is in principle no worse than the first (in 
fact, the elements are well-formed, and no triangles are needed). To serve as a basis for 
comparison, the 410 element was run with both grids, again shown in Figure 1 (curves 
labeled polar for the polar grid, and old for the regular 410 in STAGS), where the ratio 
of displacement under the load to the converged solution is plotted as a function of the 
number of nodes along a side. Two things are striking about these results: (1) that the 410 
element does so well with the polar grid, and (2) that it does so poorly with the seemingly 
harmless ‘nonstandard’ grid. This was all the more disappointing in view of the fact that 
for some coupled-field problems (such as in fluid-structure interaction applications), the 
latter grid is much preferred. One thing we noticed right away is that this problem is 
inextensional, an application where the linear bending theory in the 410 is well-suited. As 
these were all linear runs, it thus is not surprising that the 410 with the polar grid does 
well. Hence, the poor performance can be explained by studying the influence of the grid 
on the solution, and in particular, that the nodes cannot lie on the surface of the sphere 
and the elements at the same time. Apparently, equilibrium (determined at the actual 
node positions and calculated from a flat element) cannot be satisfied unless the much 
higher membrane stiffness comes into play. Thus the displacements are correspondingly 
reduced and the 410 misses the mark. 

It was at this time that we realized that our corotational theory in principle removes 
most of the rigid-body motion before the elements use the displacements [7]. Could it be 
possible to solve this same problem for a very tiny load (assuring that real nonlinearity is 
absent) using a “nonlinear” algorithm with the corotation formulation? For the polar grid, 
the solution converged in one step with an absurdly tiny error to the same linear solution. 
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But for the other grid, convergence was difficult, with several step cuts. This solution is 
labeled coro in Figure 1, and it can be seen that the results were dramatically improved. 
No matter how small we made the load, the convergence was just as difficult, and the same 
solution (except for the load proportionality factor) came out. For this reason, we began 
to look at some of the details of the corotational implementation, to see if we had taken 
all factors into account, and to see if any effects of corotation could be extended to linear 
solutions. If it were possible to correct this deficiency in a manner largely external to the 
element, another “element-independent” method could be derived to extend and enhance 
existing inexpensive elements with relatively little effort. This is the subject of the next 
section. 


2.1. The Corotational Theory: What’s Missing? 

It was apparent for some time that terms were missing in the previous formulation. To see 
this, we first summarize what steps must be undertaken when the element-independent [7] 
corotational procedure is used as part of a nonlinear iteration procedure: 


Step 1) 

Step 2) 
Step 3) 

Step 4) 


Prior to the calculation of stiffness or internal force, a new local element frame is 
calculated using the latest positions of the corner nodes. Given this new element 
frame and the latest orientations of the nodal surface triads, new deformational 
displacements are computed, including the element deformational rotations. 

Strains/stresses are calculated using the results of step 1), from which internal 
forces and (if needed) a hew tangent stiffness are computed. 

New solution increments are obtained from the linearized problem in step 2), 
using standard matrix assembly and solution procedures. The unknowns are 
incremental displacements and rotations at the specified nodes. 

The solution increments are composed with the previous solution to produce a 
new solution. The composition law for the translations is the simple addition of 
the displacement increment with the proper value from the previous iteration. 
For rotations, a new orthogonal rotation matrix is first created, and the nodal 
triad is updated by the product rule for composing rotations. 


In view of the behavior described in the previous section, we have concentrated on 
Step 2), that is, how a first and second variation of the strain energy is consistently derived 
given the element frame and deformational displacements from Step 1). All succeeding 
arguments assume that we are working within a single element, between a given iteration 
and a previous iteration. We have available to us the current solution, complete with 
deformational displacements and a new element frame. In addition, we know what the 
starting element coordinates, the starting element frame, and initial nodal surface triads 
were. Therefore, in what follows, we shall omit any reference to the assembly process and 
not worry about inter-element compatibility and assembly matters. The only reference to 
node number will be among those in a single element, a notation to be omitted when the 
meaning is obvious. Reference to the solution state will be with a single integer subscript 
in parenthesis, to be used only when necessary. 
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We now for the moment postulate that the local element strain energy is a scalar 
function of the local element deformational displacements only: 

V = (2.1) 

where U is the strain energy integrated over the element expressed as a function solely of 
quantities u d e e f defined at the nodes. For STAGS, this is at most a fourth-order polynomial 
in the displacement components. The relations 

a uf e/ = Eft) (“u, + “X, - X) - X, (2.2) 

and 

°E>(*) = E(’ fc) a S(*)E( 0 ) (2.3) 

define the deformational translations and rotations, respectively. Here, we are referring 
to state ( k ) at which a current solution is known, and to state (0) whose data we have 
saved for reference. The subscript g means that these quantities are in the global reference 
frame, and the subscript e refers to quantities defined in the element frame. X are initial 
nodal coordinates, defined relative to node 1 (see [7], where the corotational theory is 
developed in detail). The left superscript a refers to node number. There are in principle 
three translations and three rotations at each node a. 

From Eq. (2.2) one can see that for the element frame defined in [7], there are exactly 
6 N —6 independent translational quantities, where N is the number of nodes in the element. 
Since all distances are measured relative to node 1 (Eq. (2.2)), the three deformational 
translations are zero there. Similarly from the definition of the element frame in [7], it can 
be shown that the x in-plane displacement at node 4, and the lateral (w) displacement at 
node 3 are also zero. Finally, the lateral displacement at node 2 is the same as for node 4. 
Note that Eq. (2.2) is simply the statement that the local deformational displacement is 
the difference between the local element relative coordinates in deformed state k and the 
initial undeformed coordinates in the original element frame. 

Eq. (2.3) is the analogue for rotations. a D(k) is the relative deformational rotation 
matrix needed to bring the nodal triad in line with the initial triad rotated rigidly with 
the local element frame. 

The following definitions will explain the notation: 

E(fc) The local element frame for state k. Components originally defined in the 

element frame are transformed into the global frame. This is an orthogonal 
transformation. There is one and only one element frame per element per 
solution state k. 

a S(fc) The amount required to rotate the initial nodal triad at node a from the 

initial state to the current state k. There are N such rotation matrices, 
one for each node a. One a S(k) is attached to each node in the system, 
irrespective of the element connectivity. It is presumed to be uniquely defined 
and available for Step 2) above. This is also an orthogonal matrix. 
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a »(k) 


V 


The relative rotation required to bring the nodal triad (which rotates in re- 
sponse to the system demands, via the rotation composition law in Step 4)) 
in line with the same triad initially coincident but rigidly rotated with the 
element frame E(jt). This is a measure of local rotation that vanishes if the 
nodes and the element frame move as a unit. a D(*) is also an orthogonal ma- 
trix. Angular freedoms extracted from this triad [7] are used in the element 
strain computations. 

A bold lower-case Roman character refers to a vector quantity. It can be 
either an axial or a polar vector. When the distinction needs to be made, it 
will be explicitly pointed out. The dimensionality of the vector will depend 
on the context. Two exceptions to this notation are allowed: (1) the vector of 
undeformed relative nodal coordinates will be written as X (with appropriate 
subscripts), and the deformed relative nodal coordinates as x. 

A bold lower-case Roman character enclosed in a single set of brackets refers 
to a skew-symmetric spin tensor generated from a tree dimensional axial 
vector by the following correspondence principle: 


r 0 -v 3 v 


2 T 


v = 


0 -v 1 


— v 2 r 1 


0 J 


(2.4) 


Here, the superscripts refer to components , not powers. 


If the energy U( u de f) is a unique differentiable function of the deformational displacements 
including rotations extracted from D(*), then we can take its derivative. By the principal 
of virtual work, the first derivative of the strain energy with respect to the independent, 
admissible nodal freedoms will be the internal force. Unfortunately, the admissible nodal 
freedoms are not u de f and D, but global total displacements and rotations. Without 
loss of generality, we take the translations as expressed in the global frame common to all 
elements (of course, if all freedoms at a given node are uniquely defined, any computational 
coordinate system would serve just as well). For rotations, we take the current orientation 
of the nodal triad as completely specifying nodal rotation data. This triad rotates rigidly 
with the node in response to the system, updated in Step 4) above. Thus we are faced 
with finding not only the explicit derivative of U, but also the derivative of the local 
displacements as a function of global quantities, as expressed in Eq. (2.2) and (2.3). 

In the next section we shall examine the relation between the global and local dis- 
placements, and we shall reveal that for certain elements, terms presently omitted in the 
corotational theory can become important. 


2.2. The Derivative of the Element Frame 
We begin by taking the variation of Eq. (2.2): 

6 a u d e ef = (£E ( fc)) r ( a u 3 + a X g - l u,) + E(k){6 a Ug - 6 l Ug) (2-5) 
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We take note that X e and X g are constant undeformed coordinates, and that all variations 
are functions solely of the global displacements. We make one important simplifying 
assumption at the outset: that all elements to be considered are invariant to arbitrary 
translations , something certainly true for the elements in STAGS. If that is the case, it 
is permissible to arbitrarily fix the displacement to zero at node 1 without changing the 
strain energy (in other words, there are at least three eigenmodes of the stiffness, linear or 
nonlinear, that are rigid body translations with zero eigenvalues). The second term in (2.5) 
is therefore simply the transformation from global coordinates to element coordinates, 
presently performed on both the stiffness and internal force before assembly. The first 
term , however, is not now included. 

Before investigating this term, we take the variation of Eq. (2.3) to obtain 

£°D(fc) = (£E(*)) Ta S( fc )E(o) + E(’ fc )^°S(*)E(o) (2.6) 

Note that in both equations, we are faced with taking the variation of an orthogonal 
matrix. It is well known that for any orthogonal matrix T, 

[$«r] y T = 6T 
[Su T }g = (*T)X T 

where [£u>r] is a skew symmetric matrix representing the infinitesimal rotations of the 
base vectors of T. Eq. (2.5) becomes 

6° u?l = (|6w£] s E ( k)) T ( a u, + • X , - ’u t ) + E t («»u, - ( l u t ) 
which reduces to 

6*u dtf = -Bj k) \6u E }g a x g + E*>Uy (2.8) 

where we have assumed that the 1 u 9 has been suppressed because the element is invariant 
to translations. Note that we have carried out the transpose and used the fact that the 
transpose of a skew-symmetric matrix is its negative. The quantity a x g is the vector of 
deformed (undeformed plus displacements) coordinates of node a relative to node 1. Eq. 
(2.6) similarly becomes 

= -E(*)) T [£«£]y a S(*)E(o) + E(fc) [6°ws]y a S(*)E(o) (2.9) 

where similar operations to those in (2.8) were done here. Note that the skew-symmetric 
matrix \6u)E\ g is a function of the current global nodal coordinates, dependent on the 
particular definition of the element frame in [7]. [6 a ws]y is, on the other hand, an inde- 
pendent variation (in global coordinates) of the nodal triad at node a; these increments 
are the outcome of the solution process to be used in Step 4) above. 

It proves more convenient to work in element coordinates, then transforming by stan- 
dard procedures to the global frame after the local stiffness and internal force contributions 
have been calculated. The skew-symmetric tensor and its equivalent vector, for orthogonal 
matrices only, transform as follows: 

[r]y = T J( [r],Tf ( (2.10) 


r ; — Ty,r<; 
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where j and i are two different coordinate systems. Using (2.10), Eqs. (2.8) and (2.9) can 
be expressed in local element coordinates as 

6 a u d e e f = -{6u E } e a x e +6 a u e (2.11) 

and 

£ a D(fc) = ([^°W5] e - [<5u/ E ] e ) a D( fc ) ( 2 - 12 ) 

where it is understood that the components of the skew-symmetric tensors are relative to 
the element basis. °z e and “u e are the relative deformed coordinates and displacements in 
the element frame. Alternate forms for Eq. (2.11) can be derived from the correspondence 
of the product of a skew-symmetric tensor with the cross product, specifically, 

[w]r = uxr=-rxu = -[r]u/ (2.13) 

for any two three dimensional vectors r and u. Eq. (2.11) can take either of the two forms 


6 a u de f = -Sue x a x e + 

= { a X e ]6tj)E + £°U e 

where Sue is in element coordinates and where 

[ 0 - a x® 

(‘*J = 


a ^3 a ^2 


a *e 




L — z; 


a- 2 a _1 


0 J 


(2.14) 


is the spin tensor of the relative deformed coordinates for node a. Note that the spin 
tensor for node 1 vanishes since the coordinates are measured relative to that node. But 
Eq. (2.14) states that the variation in the deformational translations is modified only by 
an infinitesimal rotation of the element coordinates about node 1. Eq. (2.12), on the 
other hand, is already in the form used for the rotation composition law, something to 
be covered in much greater detail subsequently. (2.12) and (2.14) taken together state 
that the variation of the local deformational translations and rotations is modified by an 
infinitesimal rotation about node 1 with the axis of rotation in the direction of Sue- 

What does this mean physically? We return to the energy U and express its variation 
in terms of the global coordinates: 


SU = (2.15) 

du3 du' g 9 9 v ' 

where for clarity of notation we substitute u for u de f . Here, f is the internal force vec- 
tor corresponding to the deformational displacements, and J is the generalized Jacobian 
transformation from global to local coordinates embodied in Eqs. (2.12) and (2.14). The 


7 



new internal force, in terms of global coordinates (but expressed in components relative to 
the element frame) is 


%r = %r + %r\ a x e ] 


dXtlE 

d u 


“f r = % ~ % * 


dus 

d u 


= % r - *X e X % 


tr * 


dWE 

d u 


(2.16) 


where a i tr and °f r are the translational and rotational parts of the internal force, respec- 
tively. dxtiEldw is the yet undetermined variation of the spin tensor u^asa function of the 
displacements u. Thus for each node in the system, the translational force is modified by 
the cross product of the current relative position coordinates of the node with the force, 
times the variation of the base vectors of e with respect to the displacements at all the 
nodes. Each of the moments is modified by the same derivative du^/du. Taken together, 
the internal forces are transformed by the quantity 


p = I - w 


through the relation 


f=P r f 

where I denotes the identity matrix. Eq. (2.18) is fully equivalent to (2.16) when 


t duiE 

'i = 

1 du 


and 


* = 


0 


_ o x 3 


a X 3 - a X 2 


0 


°X 1 


a x 2 -“x 1 


Is 


(2.17) 

(2.18) 

(2.19) 


( 2 . 20 ) 


for a given node a. This relation is repeated for each node a, with the skew-symmetric 
portion belonging to the translations and the 3X3 identity I 3 to the rotations. 0 is a 
matrix containing a column of length six times the number of nodes N , or 6 N for the 
general case of six freedoms per node. The first column corresponds to an infinitesimal 
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rotation about the x axis, the second about the y axis, and the third about the z direction 
in the local element frame. Each of these rotations in turn depends explicitly on the nodal 
positions at the corners. Close examination of the product <j> f shows it to be equivalent 
to the relation 

N 

£(•* x “f lr + “f r ) (2.21) 

a= 1 

These are none other than the equations of moment equilibrium for the element, something 
that should be satisfied for all values of nodal displacements. It was one of the fundamen- 
tal assumptions of the corotational theory that each element will be self-equilibrated with 
respect to the current configuration. Ideally, then, these new terms, with all the compli- 
cations they bring, should vanish identically. The fact that they don’t presents both an 
obstacle and an opportunity that have motivated the continuation of this research effort. 
The derivation of this matrix 7 is the subject of the next section. 


2.3. The Variation of the Element Frame 

The infinitesimal rotation 7 of the basis vectors of E (the columns of the transformation 
from the element local coordinates to global coordinates) is determined uniquely by the 
current nodal positions for both triangular and quadrilateral elements. There is one column 
of length 6 N for each of the three directions, populated by the contributions from the 
motion of each of the corner nodes. For triangles, the element-local coordinate system has 
its origin at node 1, its y axis along the line joining nodes 1 and 3, and the z pointing out 
of the plane of the triangle. The x axis completes a right-hand system. For quadrilateral 
elements, the z axis is first defined as the cross product of the principal diagonals (lines 
joining nodes 1-3 and 2-4, respectively), the y is the projection of the 1-4 line on the plane 
with z as its normal and containing nodes 1 and 3, with the x axis completing a right-hand 
system. Details and pictures of these two cases are found in [7]. The total variation of the 
element frame with respect to the displacements is the product <fa T . This is used in Eqs. 
(2.17) and (2.18) to calculate the internal force. In the two subsections that follow, 7 will 
be derived for the special case of the triangle and quadrilateral elements, respectively. 


2.3.1. Variation of the Element Frame — Triangles 


For both triangles and quadrilaterals, only 7 need be derived as a function of the in- 
finitesimal local displacements, since the components are already in the element frame, 
and since transformation to any other frame requires only standard matrix operations us- 
ing orthogonal matrices. Recall that the unit basis vector e 2 of E lies in the direction of 
edge 1 — 3. 



( 2 . 22 ) 


For simplicity, we have omitted the reference to the element frame, as all quantities are 
defined in this frame. We recall also that these are the deformed coordinates for the 
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specified node, relative to the origin at node 1 . e 3 is defined as 


2 x X 3 X 

e 3 = —5^-. (2.23) 

Dz is the magnitude of 2 x x 3 x needed to make e3 a unit vector. Its value is, specifically 

Dz = 2 x x 3 x = 2 x 1 3 x 2 (2-24) 

|X?3 1 is equal to twice the area of the element. The right superscript refers to components 
along the local element basis vectors e*, k = 1,3. The updated coordinates “x are related 
to the undeformed coordinates by the formula 


°x = a X + a u - *u 


(2.25) 


where we must now include the displacement 1 u at node one because its variation is 
permitted. We next vary each displacement in turn to see what change in each of the three 
basis vectors is generated, in effect generating the quantity directly. For example, from 
(2.22) and (2.25) 

6*2 — *\(6 3 u z - ^u 3 ) -|- e3(6 3 u x - tf 1 u I ) (2.26) 

Taking note that the relation between [Sue]* and the variation of the element frame in 
element coordinates is (see Eqs. (2.7) and (2.10) and also Ref. [12]) 


[6o> £ ]e = E^-)tfE(i) 


we have 


and 


3 T 6 1 u 1 - 6 3 u l 

6u% = -e{ 6*2 = 3I2 


C 1 Tc 6 3 u 3 -6 l u 3 

6 u> e = e 3 6*2 = 3— 


In these two equations, we note again that the displacement components have two super- 
scripts: the left denotes the node number and the right the component number. For the 
base vectors e*, the subscript is the base vector index denoting the direction it points in 
the element system. Each column in E is one of those base vectors, and the components of 
these base vectors give the orientation of each of the three base vectors in the global system, 
since E is the transformation from element-local to global coordinates. The remainder of 
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the derivation is tedious but straightforward, yielding for 7 the relation 


0 


0 



0 0 
_i_ JLfi _ i si 

3 X 2 2 x l tl 3 r 2 

0 3 


0 

0 


0 


7 = 


0 

0 


0 

0 

1 

3 X 2 


0 0 
0 0 


1 

'**« 

0 3 


0 

0 

2 x 2 

2^13 

0 3 


0 



(2.27) 


where O3 is a three-by-three block of zeros corresponding to the rotational freedoms, verify- 
ing that 7 is a function of the nodal translations only. Translational invariance is preserved 
even though it was not included in 7 (we have assumed the elements are already invariant 
to translations): the column sums are identically zero . That is, 7 T annihilates any pure 
translation. Note also that 

^ r 7 = 7 r ^ = I 3 (2.28) 

a very important bi-orthogonality relation whose implications will be discussed later. 


2.3.2. Variation of the Element Frame — Quadrilaterals 
For quadrilaterals, we have from [7] the following: 

3 xx ( 4 x- 2 z) 

= — D, 


(2.29) 


where 

d 3 = V(V - V) - V(V - V) 
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is the normalizing factor. |Z> 3 | is equal to twice the area of the element projected on the 
plane defined by e 3 , nodes 1 and 3. We also have 


ei = 


*x x e 3 
D 2 ’ 


(2.30) 


where 


D 2 = V 

The e 2 axis required to complete a right-hand system is 

e 2 = e 3 x ei. (2.31) 

These manipulations are not unlike those for the triangles, yielding the following final 
results for 7: 

0 0 i 


h = 


h = 


7 3 = 


d 3 

D 3 


0 3 

■ 0 

0 

0 

0 

3 x‘ 
D 3 

3 z a 
D 3 


0 3 


0 

- a x» 

0 

4 z a - a x ; 

*3 

i ?3 


0 3 

■ 0 

0 

0 

0 


4 *V» 2 - 5 » a ) 

D3 


0 

0 


4 z 33 z a 




14 = 


0 

0 

4 z 3 ( 4 z a — a z 
Da D3 


1 


D3 D3 D2D3 


(2.32) 


(2.33) 


(2.34) 


(2.35) 


O3 
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where we have broken 7 into contributions belonging to each of the four nodes. This 
quantity also applies to elements with more than 4 nodes, but the contributions belonging 
to any but the nodes that define the frame are identically zero. Again, 7 does not spoil 
translational invariance, and the important bi-orthogonality relation (2.28) is obeyed (try 
that one out by direct substitution!). 

2.3.3. Variation of the Element Frame — Beams 

For beams, the definition of the element frame is more complicated. The x axis lies along 
the line joining the two beam nodes, and the z axis is the cross product of the y base 
vector (q 2 ) of a triad rigidly attached to the surface triad 1 S(^) for node 1, but initially 
coincident with the initial element frame E( 0 ): 

Q(fc) — ^(AO^qo) (2.36) 

where the initial element frame E(o) is chosen in some unique way as function of translations 
only. See Ref. [7] for details. Using the same argument as for triangles and quadrilaterals, 
the following relations hold for the local rotations 6u>e< 

6u>g = —e^Sez 

Su% = -el6e x (2.37) 

E = ® 2*^®1 

The variation of the base vector of the element frame ei is a function only of displacements, 
and is derived in exactly the same way as with triangles. The result is 

6ei = y [e 2 (£ 2 u 2 - 6 l u 2 ) + e 3 (6 2 u 3 - 6 1 u 3 )] ^ 

/ = V 

where u are displacements with components labeled by right superscripts and the node 
designation by left superscripts. / is the current length of the beam element, which by defi- 
nition lies along the x axis. The variation of e 3 is more complicated, but is straightforward, 
made much simpler by using the relation 

e 3 6e 3 = 0 

After tedious calculations, the result is 

1 2 

6e 3 = - «‘« 3 )] - V - <‘u 3 )} (2.39) 

92 * * 

including all terms to first order. Note that all quantities are defined in the tltmtnt frame, 
as was the case for triangles and quadrilaterals. The expression for 6w% is derived by 
substituting (2.39) into (2.37) to yield 

<«I = 4w-t(‘ v - <, “ 3 )! ( 2 - 40 > 

92 * 
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6Q e is related to the variation of the surface triad by the relation 


6Qe = [ u s]eQ e 


(2.41) 


where we have dropped the system state information and inserted the subscript e to empha- 
size the use of element coordinates. Combining (2.40), (2.41), (2.38), and using the same 
notation for 7 as we had for triangles and quadrilaterals, we come up with the following 
final result: 

■ 0 0 0 ' 

0 0 -} 


1 0 0 


-0 0 0 

0 0 0 


(2.42) 


0 
0 

_ 1 

l 

03 

4 

Note that if for some reason the beam is so deformed that q 2 lies along the z axis, our 
algorithm switches to the choice of the q 3 as our e3 base vector, in which case everything 
here is derived as before except that 0 is set to zero. We leave that one as an exercise for 
the reader! 

To calculate 7, we need q 2 in element coordinates. To do this, we use the relation 
(2.10) (which applies to orthogonal matrices as well as spin tensors) and extract the second 
column: 

q 2 = E(’ fc ) 1 S(fc)E( 0 )e 2 (2-43) 


0 

0 

l 

where 


0 


1 

1 


0 


It is gratifying that the bi-orthogonality relations still hold between <f> and 7. The 
principal difference here is that the variation of the element frame is no longer independent 
of the rotations, with contributions dependent on the variation of the surface triad. 
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2.4. The Projection Operator P 

Up to this point we have demonstrated that the matrix represents an infinitesimal 
rotation about node 1 in three orthogonal directions corresponding to the base vectors 
of the element frame E. We have also shown that the inner product of <f> with the local 
internal forces reproduces the violation of rotational equilibrium in the three directions (Eq. 
2.21). What remains is to show that the effect of infinitesimal rotations on the internal 
force modified by the product 70 r vanishes. If the element with modified internal forces 
is to satisfy equilibrium in the updated configuration, it must be invariant to infinitesimal 
rotations. Invariance to rotation is completely equivalent to the satisfaction of equilibrium, 
or the requirement that the element taken as a free body will be in rotational balance and 
will not move. What we must show is that 

0 T P T f=O (2.44) 

where, 

p = i - <h T 

But given the bi-orthogonality relation (2.28), Eq. (2.44) holds. In fact, if a general 
displacement d is decomposed into 


d = &dtf + 0® (2*45) 

where d <* e / is the nonrotational part of the displacement, and a is a vector of weighting 
coefficients for the rotational part of the displacement, we have 

Pd = (I-^, T )(d d , / + «a) 

= <W + d[is - b T d)]« - <h T &i.s 

The only way the contribution from the rigid-body part is always zero is if 


7 T 0 = Is 


which is identical to the bi-orthogonality relation (2.28) that we have shown to hold for 
7 and 0. Thus the modified internal forces are strictly invariant to infinitesimal rotations 
(and translations for the elements in STAGS). In addition, if the element nearly obeys 
equilibrium in the first place, the quantity <frf will be small, and the internal force and 
stiffness for the deformational part of the displacements will be only slightly altered. This 
translates into improved performance for the warping problem, because it is the rotations 
about an axis perpendicular to the element plane, which can be much larger than de- 
formational displacements across the area of the element, that spoils performance. In a 
sense, this is the analogue of the corotational theory in the linear limit, because if local 
equilibrium fails in that limit, the influence of P does not vanish. 

Before reviewing the effect that P has had on various elements, we first mention that 
P is a projection operator. Projection operators have no inverse (this P is of rank N - 3), 
and satisfy the relation 

P 2 = P 
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Projection operators are annihilators: they project out or reduce the vector space, in 
this case, removing the three rigid rotations. What is powerful about this particular 
projection operator is that it was consistently derived from the corotational theory that 
has already proven itself. The vectors <f> and 7 are well conditioned (unless the area 
of the element is zero) and dimensionally consistent, with the translational parts of <j> 
having the units of length, and 7 its natural inverse with units of one over length. No 
arbitrary orthogonalization was necessary, as orthogonality just “fell out” as a result of 
taking consistent derivatives. This suggests that the method will be applicable no matter 
what element frame is chosen, such as the more complicated element frame for beams [7]. 

Finally, we mention that while the internal force is modified by the relation in Eq. 
(2.18), the linear stiffness is modified by the usual equivalence transformation 

K = P t KP (2.47) 

The derivation of (2.47) is straightforward, and its implementation follows standard pro- 
cedures already in place. 


2.5. Element Performance — Results 

In the beginning, the projector was added to the 410 element and tested with planar and 
cylinder configurations that do not have warping. No difference whatever was seen between 
the old and new 410 for linear analyses, indicating that the 410 obeys equilibrium in the 
linear limit regardless of element distortions within a plane. The result reinforces the belief 
that, for a majority of cases, the current implementation of the 410 suffices. 

Figure 1 shows what happens when the SH410 element is run with the projection P in 
place (linear analysis) for a case where warping does exist. Note the improved performance 
of the 410 even over the nonlinear results discussed earlier. As expected, this time when 
the projection operator was included, the nonlinear solution required only a single iteration 
(with an absurdly small error) that reproduced the linear solution almost exactly, even at 
the unit loads used in this analysis. 

Figures 2 and 3 summarize the effect of P on various elements. Included are the 
STAGS 410, the four node ANS [6] (4.ANS), an experimental hybrid element [4] (HYB), 
and the nine node ANS curved element [5] shown for reference (it is curved and is insensitive 
to warping). Fig. 2 is the same pinched hemisphere problem as for Fig. 1. The 4.ANS is 
K. C. Park’s element that has been shown to be sensitive to warping. Clearly the projector 
improves performance significantly. Even the HYB element is helped, which means local 
equilibrium is not truly satisfied. 

Fig. 3 is the pinched cylinder problem presented in [10], where this time distortion 
angle is the parameter. The grids chosen were different for each of the elements, so that 
for zero distortion, one begins with a converged solution. The sensitivity to distortion (a 
typical grid shown in the figure) is measured by computing the solution as a function of 
the orientation of the central element (see figure). At about 40°, the area of the adjacent 
elements is nearly zero, presenting a “worst case” . Warping is important in this problem 
because when the grid lines no longer follow the generators of the cylinder, the element 
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nodes cannot both lie on the cylinder and on the flat element segments. The results for 
this problem are just as dramatic as for the pinched hemisphere. All flat elements are 
improved, with the STAGS 410 element one of the best performers. What a contrast to 
the unaltered case! The old 410 essentially locked and was completely useless for even the 
most modest distortion angles. We want to mention that for the case of r/h = 50, all 
elements reproduced the converged solution no matter what the distortion. These results 
are not shown here. For a look at how they performed before projection see (10] . 
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3. The Problem of Slow Convergence 

Soon after the corotational theory was implemented in STAGS, we recognized that con- 
vergence to nonlinear solution states had deteriorated somewhat, depending on the nature 
of the case investigated. True Newton cases no longer acted like true Newton, in that 
the quadratic, or superlinear convergence was lost. Clearly the second variation was not 
completely consistent with the first variation of the strain energy; this problem appeared 
whenever corotation was turned on, regardless of the type of analysis, and regardless of 
whether the external loads were conservative or not. A clue that some degree of incompat- 
ibility that might exist in the tangent stiffness was clear even for the Mak arch [3] chosen 
in [7] as a benchmark for the corotational implementation. For the SH410 element, the 
sign of the determinant of the tangent stiffness did not change at the instant the load 
crossed the limit point, as of course it must. An associated problem was the sensitivity of 
the corotational analysis to the frequency of element frame updates. The solution came 
out differently when the frame was updated every iteration as opposed to each converged 
load step as was originally proposed in [7]. This was disturbing because the 41X family of 
elements was advertized and tested for moderate rotations at least up to 10 or 15°. Ideally, 
one should be able to “turn off” corotation for a while, do some load steps in the neigh- 
borhood of a converged state using as a base the element frame calculated there, and then 
repeat the process. One could realize possible improvements in efficiency, and certainly 
have better confidence in the method if the answers were not sensitive to such choices. 
The fact that, as corotation is now configured, the answers are sensitive to the frequency 
the element frame updates has been part of the reason for conducting the research that 
has led to the results to be covered in this section. 

Section 3 is divided into three main subsections as follows: 

1) In the first subsection, we shall demonstrate nonsymmetry in a numerically- 
calculated tangent stiffness (using only first variations), even for a supposedly con- 
verged equilibrium state. We shall discuss what consequences this nonsymmetry 
has on both the underlying assumptions of the corotational theory and solutions 
calculated using the theory. 

2) In the second subsection, we shall introduce three simple sample problems to explain 
what is going on. Each problem has an analogue in the real finite element world, and 
together they cover the features we are interested in, including the nonsymmetry 
questions in 1). 

3) In the last subsection, the effect of the second derivative of the relations (2.2) and 
(2.3) that define the transformation from the global, total displacements to the local, 
deformational displacements will be discussed. A way to account for the greater 
contribution of these derivatives based solely on global invariance properties will be 
derived. 

By combining the what is learned from these three areas, we will demonstrate the recovery 
of true Newton convergence for a real STAGS 410 analysis for a problem with conservative 
loads. 
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3.1. The Nonsymmetric Tangent Stiffness 

After we had derived and implemented the projection operator P for the SH410 element, 
the next question was what would be the effect on a nonlinear analysis. We first altered 
the material (as opposed to geometric) stiffness by the relation in Eq. (2.47), and tried this 
on some sample test cases. We were gratified that the results did not change much, except 
that convergence was no better than without P. We discovered that P did, in fact, modify 
the nonlinear results slightly even for planar problems, indicating that the geometrically 
nonlinear element no longer obeys equilibrium relative to the deformed state. 

The next logical step was to get a closer look at the tangent stiffness for the element in 
some nonlinear configuration that contained nontrivial rotations. One excellent means of 
doing this is to do as many first variations as there are freedoms, perturbing each freedom 
in turn by a “tiny” amount, and taking the ratio of the difference from some reference 
state to the increment chosen. As this is the definition of the partial derivative, it should 
converge for some value of the increment. It is to be noted that by perturbation, we mean 
traversing all the steps in the corotational theory from Step 1) through 4), omitting only 
3) because the increment is a given and not solved for. This means that when a rotational 
freedom is incremented by a small amount, a full rotation matrix corresponding to this 
increment is created and composed with the previous nodal triad to produce a new one. 
Rotations are combined by the product rule as follows: 

a S(fc+i) = A a S (fc+1(fc ) a S( fc ) (3.1) 

where we are using the proper notation for the nodal surface triad. The quantity 
A°S(fc +1) k) is the rotation matrix, presumably very close to the identity, derived from 
the perturbation of a given rotational freedom. This rule is in contrast to the ordinary 
addition rule used for translational freedoms. 

When this process is carried out (including the calculation of new element frames, 
etc.) the resulting stiffness turns out to be nonsymmetric and different in many respects 
from that calculated by standard means. This poses a problem, because it is well known 
that for a conservative system loaded only by point forces in a fixed direction, the stiffness 
matrix must be symmetric if the linearization is carried out consistently. At this point we 
decided to backtrack to some simple problems that embodied Steps l) through 4) in the 
full corotational theory, but at the same time could be fully understood. 


3.2. Sample Problems — Introduction 

We will begin with a problem defined by unit position vectors only, a problem whose 
linearization is straightforward, easily understood, and important in its own right. We 
shall follow this with a problem that is instead a function of a pseudovector (angles), 
and show what has to happen when the product rule is used to update the system state. 
Finally, we shall investigate what difference the choice of the pseudovector makes (if any), 
and what this means for a finite element analysis. The third and the second problem will 
be the same except for the normalization of the pseudovector. 
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3.2.1. The Function of a Unit Vector 


Our first sample problem has a potential that depends on unit vectors only, including 
the external force. This problem belongs to a general class of potential functions where 
rotations are replaced by shell normals at the nodes, so that exactly 2 rotational freedoms 
exist per site. Most C° shell elements belong to this class, where there are five freedoms 
per node, with three translational freedoms and two rotations. Potentials defined by these 
normal vectors are well-defined, conservative, and very straightforward to linearize. 


3.2. 1.1. The Definition of the Potential U 
The sample potential chosen is 

1 I 3 

U{fi,r 2l r 3 ) = -{kiel + k 2 el + k i2 e 2 12 ) + (3.2) 

6 t=i 

where the index i is counted cyclicly from 1 to 3. For example, if * = 2, then * + 1 = 3 and 
i + 2 = 1. The subscripts label three mutually perpendicular vectors ? t that rotate rigidly 
as a group. The roman superscripts denote components of ? . The e’s are “strains” defined 
by 

f i = f? = e 3 • r, 

62 = f 2 (3.3) 



where we have introduced the notation e* to denote the base vector in the fixed coordinate 
system. This problem is like a little “beam” segment tied to springs (Fig. 4). The first 
three terms in U simulate the bending energy of the segment. The second part of the 
potential defined by the sum simulates external moments to which they reduce for small 
angles, k and m are constants that define the magnitudes of the “internal energy” and 
“external forces,” respectively. 

3.2. 1.2. The Linearization of the Potential 

Eqs. (3.2) and (3.3) constitute a conservative, nonsingular system with three freedoms, 
that is, the three quantities that define the rotation operator that brings this rigid triad 
{r 1 ,r 2 ,r 3 } in into a position of minimum energy. These three freedoms are defined implic- 
itly through the nodal triad, just as in STAGS and other structural analysis codes. There 
is no ambiguity if we begin by selecting some initial position for triad ri,r 2 ,r 3 , and by 
selecting the identity as the initial rotation operator. After the initial state is specified, 
we simulate Steps 1) through 4) (Section 2). During this entire sequence, the freedoms 
appear explicitly only as increments during the solution and update stages (Steps 3) and 

4 )). 

Steps 1) and 4) can be combined after the initial iteration because after the update, we 
have available the “local” rotations from the triad {ri,r 2 ,r 3 } that completely determines 
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the energy and its derivatives. For Step 2), we must linearize the energy with respect to 
infinitesimal rotation increments. Fortunately, for vectors, we know the relation between 
an infinitesimal rotation and the variation of the vector: 


6rk = Suj x rk 

for any small rotation 6u. In what follows, we will need to use the relations 

uvxw=uxvw 
U X (v x w) = v(u • w) - w(u • v) 
u X V = — V X u 


(3.4) 


(3.5) 


from vector analysis. Using (3.4) and (3.5) and the notation e* to denote unit fixed base 
vectors in the global system, we have 


/Cl o£l 9 

6U = [Aifjfi X e 3 + fc2 c 2^2 X ©3 + — “ (fi X ©2 — ?2 X ej) 


! 3 

- X ®t+ 1 -r.'+i X e 1+2 )] ■ 6 u 


(3.6) 


1=1 


where the same cyclic summation over indices is used. 

The second derivative is straightforward, but tedious, with liberal use of relations 
(3.5). This tangent stiffness separates naturally into two parts, the first being a “ material ” 
stiffness wherein the variation takes place again as a function of the unit vectors, to be 
then expressed as a function of the incremental rotations 6u. Specifically, 


t ,„ t-T &U (dU 

i U = <r - » w/ ri + 6li Wi 


(3.7) 


where the first term is the “material” stiffness. The second term is a “ geometric ” part, 
defined by the second variation of the unit vectors r* as a function of the incremental 
rotations. The material stiffness is 


K mat = *i(*i x ®s)(*i x e 3 ) r + M * 2 x e 3 )(f 2 x e 3 )' 

+ 7 (*i 2 (ri x e 2 - r 2 x ei)(fi x e 2 - r 2 x e x ) T 
4 


(3.8) 


Note that the “external loads” do not appear in the “material” stiffness, and that the 
material stiffness is strictly symmetric, regardless of the state of the system, equilibrium 
or not. The “geometric” stiffness is 


1 3 

K geom = - y^m,[f, + 2e?+i - et+ie^f2 + * e »+2 - ft+2e«+i)J 


i=i 


+ - Ifi • e 3 ) + - I? 2 • e 3 ) 

^12fl2 r * _T 


(3.9) 


+ 


fief - fef + I(r 2 • ei - r i • e 2 )] 
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where cyclic index notation applies, and where we have kept the inner products r^ey = r J k 
to emphasize the structure of the geometric stiffness K seom . 

The most obvious property of the geometric stiffness is that it is, in general, nonsym- 
metric. However, it can be shown by direct substitution that if (3.6) holds, the stiffness 
is always symmetric. That is, for any equilibrium configuration, the stiffness must be 
symmetric for conservative systems. Note that by conservative, we mean definable by a 
potential. The first and second variations of the “external” loads must also be consistently 
derived, which of course includes the contribution of the so-called “load stiffness” to Eq. 
(3.9). The second obvious property is the repetition of the sequence 

r k ej - If k • ey (3.10) 

for various indices k and j. It turns out that there is a generalization of (3.8) and (3.9) 
that makes use of this pattern. We shall show this subsequently. 

3. 2. 1.3. The Solution of the First Sample Problem 

Small, in-core utilities that simulate Steps 2), 3) and 4) were put together into a true 
Newton- Rap hson solver to be used for arbitrary combinations of the constants A:, and m,-, 
and arbitrary starting configurations of the triad {r!,^,^}. The surface triad update 
routine from the STAGS corotational software (Subroutine UPSTR) was altered to create 
a rotation operator from the relation 

S(ws) = exp([u;s]) (3.11) 

This is the general relation connecting any skew-symmetric matrix [ws] to the rotation 
operator S. Eq. (3.11) has been discussed at length by ARGYRIS [2], SIMO [8,9] and 
many others. The modifications to UPSTR are minor, involving only the renormalization 
of the pseudovector defined in [7] to conform to (3.11). In all other respects, the solution 
sequence was just like any ordinary Newton-Raphson iteration loop, with the definition 
of the internal force by Eq. (3.6), and the tangent stiffness by sum of K ma t and K geom 
(Eqs. (3.8) and (3.9)). A new factorization was carried out for each iteration, solution 
increments were solved for, and Eq. (3.11) was invoked to create an incremental rotation. 
This incremental rotation was applied to the latest configuration of the triad {fi,? 2 ,r 3 } 
to create the next configuration. The error norm of the residual (Eq. (3.6)) was then cal- 
culated for output. This cycle was repeated as many times as desired, and the convergence 
was manually monitored as the solution progressed. We have included the option to leave 
the stiffness alone, or to symmetrize it before factorization was put in place to see what 
effect this would have on convergence. 

The first runs were made using the ordinary, nonsymmetric stiffness. As many arbi- 
trary combinations of the input were tried as time permitted, with the outcome always 
exactly the same: after a few starting iterations (that number depended somewhat on the 
initial conditions), the error diminished at least quadratically like true Newton should, with 
the error norm at step i being roughly the square of the previous norm. This behavior 
shows that the linearization in Eqs. (3.6) through (3.9) is completely consistent with the 
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update formula (3.11), and that (3.4) indeed defines properly the variation of a unit vector 
in three-dimensional space. The next runs were identical, except that the stiffness was 
symmetrized first. The outcome to this set of runs was almost identical! It was even the 
same when the transpose of the stiffness was used, implying that any scalar times the 
skew-symmetric part of the stiffness could be added to the symmetrized stiffness without 
substantial change to the solution convergence. Bahram Nour-Omid has explained this 
behavior with a very simple proof, shown in the next section. 

3. 2. 1.4. Convergence with the Symmetrized Stiffness 

Consider the semi-Newton algorithm given as follows for the solution of the equation 

f(x) = 0 (3.12) 


where f is a vector function of vector x. 

Given an initial guess x(t): 

1) K (0 = U 

X=X <«> 

2) Solve K; 0 Ax (0 = — f(x (0 ) (3.13) 

3 ) 

where we define the quantities 

K* 

K tt 


*(.■+!) = x (0 + Ax (») 


^(K + K T ) 

hK-K T ) 


(3.14) 


Theorem 

The rate of convergence of the above algorithm defined by steps 1) through 3) in (3.13) as 
constituting an iteration cycle is quadratic if 

l|Kf„|| < a||f(x (l) )|| (3.15) 

where a is a constant. 

Proof 

The Taylor series expansion at X(*) is 

f( x (»')) = f ( x (*'-i)) + X(i_i)Ax(,-_i) + r(Ax(,_j)) (3.16) 

where r is the remainder which depends on AX( t _ X ) and 

l|r(Ax (j _ 1) )||</J||Ax (j _ I ,|| 2 (3.17) 
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with 0 some other constant. Writing K( t _!) in terms of its symmetric and skew-symmetric 
parts, we have 

f(x (t) ) = f(x (t _ 1} ) + Ax (t _i) + K“ t _ 1} Ax (t _ 1) + r(Ax ( i_,)) 

By using the second of Eq. (3.13) for the i — 1’st iteration, we have 

f(x (l) ) = K® t _ 1) Ax ( ,_i) + r(Ax ( ,_ 1) ) (3.18) 

Then, 

Ax (t) = -(K(, ) ) -1 {K®,_ 1) Ax(,_i) + r(Ax ( t_i))} (3.19) 

and taking norms 

IIAxwll < ll(KJ 0 )- , ||{||Kf j . 1) || ||AX(j_ 1) || + l|r(Ax ( i_,))||> (3.20) 

But 

l|K?.- l) || < o|Wx (< _ 1) || < o||(KJ < _ 1) || ||Ax (( . 1) || (3.21) 

Then, 

l|A*(oll < ll(K(.))' 1 ||{o|l(K; j . 1) )- , || ||Ax ( ,. 1) || s +^||Ax (i . 1) || 2 } (3.22) 

or 

I|X(< +1 )-X (< )|| <7i||x (t ) -x^ijll 2 (3.23) 

where 

H ~ I l(K a )~ 1 1 |(ot|J(K*,-_ 1) )— 1 |j + (3.24) 

Q.E.D 

Thus if the antisymmetric part of the stiffness vanishes as the same rate as the residual, 
convergence is quadratic whether the symmetrized or unsymmetrized stiffness is used. But 
this is just our situation: the nonsymmetric part of the stiffness is proportional to the 
residual, something to be proved for the general function of a unit vector, in the next 
section. 

3.2. 1.5. The General Function of the Unit Vector 

We now generalize what we have found for the special, but reasonable potential defined by 
(3.2). Take any function of a set of unit vectors a f, and for the time being assume that 
each is independent of the other (two independent freedoms defined per vector): 

U = f7( 1 r, 2 f,---, N f) (3.25) 

defining some twice-differentiable function U as a scalar dependent only on the components 
of the a r for “nodes” a = 1, N. The first variation, using (3.4) and relations (3.5), for each 
“ node ” a is the following: 

a f-£ a u> = °r x °f-£ a u; (3.26) 
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where 


(3.27) 


is a vector of length three belonging to “node” a. Eq. (3.26) has an obvious interpretation: 
“f is the “moment” constructed from the cross product of the unit vector with the “force” 
a f at its tip. The condition of equilibrium is for these internal moments to balance any 
external moments. The residual from an assemblage of moments so calculated must also 
vanish. The problem defined by (3.2) is such an assemblage. We in that case had three 
vectors driven by a single rotation. The first variation of (3.2) could have been derived by 
first taking the derivative with respect to the components of the three vectors as if they 
were independent quantities, then using (3.26) to relate them to three separate rotations. 
Since these three rotations are actually the same, the parts belonging to a given component 
of these rotations are added, yielding the same results as those computed directly ( Eqs. 
(3.6), (3.7) and (3.9)) for the first and second variations. Of course this step is completely 
analogous to the assembly step in finite element analysis (part of Step 3), Section 2). 

For the second variation, it is straightforward to derive the following: 

K mot = Z r KZ (3.28) 

where Z is block diagonal, with exactly N blocks, or one for each independent rotation. 
Each block contains skew-symmetric 3x3 matrix 

a Z = -[“*] (3.29) 


The second part of the stiffness is 


a K se0 m = °T a i r ~ IaC? • “f) 


(3.30) 


Note the similarity of (3.30) to (3.10); the fact that (3.9) is a special case of (3.30) shows 
clearly as a repetition of the structure in (3.30), with N of these °K 9eom blocks placed 
along the diagonal. (3.28) and (3.29) are derived by considering a single-node system and 
generalizing. The total stiffness is the sum of the “material” and “geometric” contributions. 
In (3.28), K is the second variation of the potential with respect to the components of the 
unit vectors, or 


K = 


a 2 # 

d a rd 6 r 


(3.31) 


The matrix in (3.31) could, in principle, be full. 

Some general comments are in order. First, each submatrix is of rank 2. This is true 
because 

[ fl r]°r = a r x °r = 0 (3.32) 


and because 

°K,eom a r = 0 


(3.33) 


Thus, a set of rotations each in the direction of its unit vectors cannot change their position, 
hence cannot change the energy or internal force. 


25 


Second, the tangent stiffness, is, in general, nonsymmetric. However, when the mo- 
ments are zero, then 


a r x a f = 0 


so that 


a f = c a r 


where c is some constant. Clearly then the stiffness is symmetric, confirming the proposi- 
tion that the stiffness must be symmetric whenever a conservative system is in equilibrium. 
Direct calculation of the off-diagonal skew-symmetric portion of the geometric stiffness 
confirms that skew-symmetric portion of the stiffness vanishes according to the following 
general relation 

•s;„ m = (3.34) 

something proved by ARGYRIS in [l] for the general case. Hence, for any system that 
is a function of unit vectors, we can use (3.26) for the first variation, (3.28) and the 
symmetrized version of (3.30) to complete the linearization with the second variation. 

For those functions that depend also on position vectors, we can easily generalize what 
we have. The first variation for the translational part does not require the modifications 
derived here; for the rotational part, the internal force is modified by the relation 


a f= a fx“f (3.26) 

for each node a. For the second variation, °Z is generalized to contain the identity wherever 
translational freedoms occur, with the remainder of the definition the same, as follows: 


Z = 


i 3 0 

0 °z 


(3.35) 


The results from Eq. (3.30) are added only to each rotational block of stiffness; only 
freedoms belonging to single node are coupled by (3.30). 


3.2. 1.6. Implications for General Finite Element Analysis 

Since there is, in principal, no restriction on the form of the potential defined by (3.25), 
what was developed in the last subsections applies also to C° finite elements that can be 
expressed as a function of unit shell normals at the nodes. Referring to Steps 1) through 
4) in Section 2, we note that no modifications to the sequence of calculations are necessary, 
provided that 
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1) The unit vector positions are calculated using the deformational rotation operators 
directly, rather than first converting to pseudovector angles, as is done with the 
corotational implementation in STAGS. The modification to the software to meet 
this objective is trivial. 

2) The linearization defined by (3.26) through (3.31) should be carried out. In most 
situations, it means only computing the symmetrized version of (3.30) and adding 
this 3x3 contribution to each rotational block. The additional computational load 
is very small compared with the existing implementation of relations like (3.28). 

3) The update formula for rotations must be modified to conform with Eq. (3.11). 
This has been done, and proves to be simple, stable, and efficient. 

3.2.2. The Function of a Pseudovector 

Although the results from Section 3.2.1 can be carried intact to many C° elements whose 
rotational freedoms can be described by unit vectors attached to each node, there are 
many other elements inherited from moderate rotation theory that use angles directly in 
the element interpolation. For the SH410 family in STAGS, there are three independent 
rotations at each node. For the ordinary addition rule for small rotation composition, 
no inconsistency exists (and true Newton convergence is observed), but the elements are 
restricted at most to moderate rotations. When the corotational theory was implemented, 
the correct objective method of updating rotations was invoked as 

a S(fc+-i) = A a S( fc+ | ) jfe) a S( fc ) (3.1) 

where we have repeated Eq. (3.1) for emphasis. With the introduction of (3.1), the 
consistency is lost, and so is the efficiency of the convergence to the solution. We now 
investigate what is really going on during a Newton iteration, and introduce a sample 
problem to check it out. 

3.2.2.I. The Linearization Process 

In all iteration algorithms, the goal is to obtain solution increments that bring a current 
configuration into a new configuration that satisfies equilibrium requirements, and to do 
this in as few iterations as possible. Our measure of satisfaction of equilibrium is the 
residual , or some vector function f(x(<)) that vanishes when we are done. The sequence 
of events is stated in Eqs. (3.13). Note that when convergence is taking place, f is very 
small. Let us assume that near equilibrium we have 

f(i) = rf(0 (3-36) 

where f is the directional derivative of the internal force with respect to the parameter 
e, with £ some tiny constant that vanishes when equilibrium is satisfied. We also have 
shortened the notation, so that the system state information subscripted to fy) means of 
course, f(X(,)). We wish to solve 

KAX(,-) = — ef(t) (3.37) 
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for some small value of e. We see, however, that on the one hand, the update of the system 
rotations follows (3.1) instead of (3.13), so we have, using (3.11), 


S«+i) = exp(e[ws](,-))S(f) (3.38) 

where }u>s] is also intrepreted as the directional derivative of the rotational part of x 
with respect to the parameter e and is the outcome of the solution of (3.37). Clearly 
the increments depend on e because of (3.37), vanishing as the residual vanishes. On the 
other hand, if the element strain energy depends instead on angles instead of a D( t ) (see 
Section 2), we must connect this variation with the what happens in (3.38) as e vanishes. 
Suppose we express °D( t ) as an exponential of the pseudovector that we use in the element 
calculations: 

a D (t -) = exp([(? £>](,)) (3.39) 

Now, presuming the element frame remains fixed, the variation in the relative rotation 
°D(») is 

*‘ D <0 = [«w s ] (i) a D (0 
= t[*s]*D (0 

which is from ( 2 . 12 ). In the limit of small e, the difference between two surface triads 
computed by (3.38) is identical to (3.40). Taking the variation of (3.39) and comparing 
this with (3.40) we have 

6{exp([tfo]( t ))} = [£«s](»)D(i) (3.41) 

Noting the equivalence of the solution increments through (3.38), we wish to solve for 
[ 6 u> 5 ](,-), as follows: 

[*«s] = ^{exp([tfD](i))}Dj) (3.42) 

What we need then is the derivative of the exponential of a skew-symmetric matrix times 
its transpose, or 

(3.43) 

for any skew-symmetric [#] dependent on some parameter e. If we have this connection, 
we have consistency with the update procedure in (3.38). 

This missing link has been provided us by SIMO in [9], which, when translated to the 
normalization in [7], is described by the relation 


r-n [*1 + KM - 

* exp H #| ) 


where the normalization in [7] is 


0 = 


tan(| 0 |/ 2)0 

1 * 1/2 


(3.44) 


(3.45) 
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and where in both (3.44) and (3.45) the dots over the symbols represent derivatives with 
respect to the parameter e. The axial vector equivalent of (3.45) is 


o+f$xd 


(346) 


where we have taken note of (3.42), and assumed that Su> is a function of c. We have 
dropped system state indices for the moment. Eq. (3.46) is the link we are looking for 
to complete the linearization for elements like the SH410. We next must invert (3.46) to 

express 0 as a function of the subject of the next subsection. 


3. 2.2.2. The Inversion of the Exponential Derivative 

It is tedious, but straightforward to invert (3.46). First, we express the inverse as 


$ = au + (30 x u + ^(O • u)0 (3.47) 

because the vectors u,$xu, and 0 certainly span the 3x3 space described by (3.46). If 
we substitute (3.47) into (3.46), we must get the identity. This is true only for 



(3.48) 


To get to this point, we had to make repeated and careful use of relations (3.5). The final 
formula we want is 

j = u>-i*xw+i(*.w)* (3.49) 


3.2.2.3. The Linearization of a Function of 0 

We are now in the position to linearize a potential that is a function of 0. Suppose we 
have such a potential defined by 


U = U( l 0, 2 0, (3.50) 

with the indices ranging over “nodes” 1 to N , exactly like the definition of the potential 
in (3.25). Letting the new “internal forces” be defined by the relation 


a f = 


dU 

d a e 


(3.51) 
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we have, using (3.49), the relations 


S a u Ta {= 6 a 9 Ta f 

= S a u T { a f + ^ a e x °f + a 9 • °f ) a 0) 

Thus the internal forces derived from angles are modified by 


(3.52) 


a f = a f + - a $ x °f + ~{ a 0 ■ a f) a 9 (3.53) 

2 4 

Great care must be taken since the order of the operations used to derive (3.52) is impor- 
tant. Close examination will show that the transpose of (3.49) is used for the consistent 
linearization. 

The second derivative is obtained by differentiating (3.53) and using (3.49) again. 
This is complicated but straightforward to yield 


K mat = Z'KZ 

a K geom = -|[ a f] a Z + j{( a *. a f) a Z + °* ( a f Ta Z)} 

a 4 


where each 3x3 block of a Z is defined by 


(3.54) 


“Z = I 3 - i(“}| + -“S ‘t T (3.55) 

2 4 

Of course, (3.55) is just the matrix equivalent of (3.49), as can be seen from Eq. (2.13). 

3.2.2.4. The Question of Symmetry 

Note the similarity to (3.28) through (3.31) that define the linearization of a function 
of a unit vector. This shows that there is a strong link between functions of a unit 
vector and functions of the pseudovector 9. The “geometric” stiffness is again, in general, 
nonsymmetric, just as for the function of the unit vector. It can be shown by direct 
substitution and comparison with (3.53) that once again the stiffness matrix becomes 
symmetric according to the relation (3.34). Thus all we have said about questions of 
convergence for conservative systems also applies here. 

3. 2.2. 5. The Test with a Trial Potential 

In order to test what we have come up with, we went through all the steps described by 
Eqs. (3.53) through (3.55) with the simple quadratic potential 

U(9) = \ (9 - 9 0 ) t K (9 - 9 0 ) (3.56) 

a 

where K is a 3 x 3, symmetric, positive-definite matrix input by the user, and where 9 0 is 
some constant starting value of the pseudovector completely equivalent to some “applied” 
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loads. Except for the relations stated in this section, the trial software was the same as 
for the function of the unit vector. We included the option to symmetrize the stiffness, or 
to use its transpose. In contrast to the previous program, we could start with the nodal 
triad S(o) in any conceivable orientation, enabling us to have the most general set of initial 
conditions. 

The results were gratifying. No matter what the initial state, and no matter what the 
input stiffness was, after a very few iterations, the convergence was quadratic, showing once 
again the complete consistency of the linearization. Any tinkering with relations (3.53) 
through (3.55), such as omitting some terms that may be considered small like those 
quadratic in 9 completely destroyed the quadratic convergence to the solution. However, 
the symmetrized matrix produced the identical favorable convergence behavior, as did its 
transpose. The solution was insensitive to the initial state except when the pseudovector 
was computed for angles near 180°, where it is known to be singular. This is of no 
consequence to ordinary finite element analysis, because, of course, the angles we are 
speaking of are the local, small, relative rotations for which only moderate rotation theory 
need apply. 

It is therefore clear that this linearization is complete, and that what was said about 
functions of a unit vector carries over intact to functions of a pseudovector. In particular, 
the symmetrized version of the stiffness works for any conservative system. The band 
structure of the matrix Z is described by Eq. (3.35). Only freedoms connected to a single 
node are coupled by K geom . The basic matrix operations are identical for the two cases. 
The only real fundamental difference is that here, for the pseudovector 9, we can have three 
independent freedoms per node, as the submatrices described here can, in general, be rank 
3. This is also gratifying, because elements using the pseudovector such as the SH410 can 
have three independent rotations with finite stiffness. 


3.2.3. Functions of Angles 

Up to this point, we have assumed the element rotations could be described by 9, which 
has the tangent normalization described in [7] and by (3.45). What happens if we use 
real angles? By real angles, we mean that if one rotates about a fixed axis by a certain 
angle, the magnitude of the pseudovector will agree with this angle. We shall call this 
pseudovector 9; it is this quantity that appears in Eq. (3.3d) for the exponential of the 
skew-symmetric matrix. We discuss the linearization next. 


3.2.3. 1. Linearization of a Function of Angles 
The equivalent of (3.46) for 9 is also taken from [9]: 


u — 


sin 1 9 

jij 


9 + 




9_ 1 

9 1 + 2 


( sin(l#|/2) 

V 1*1/2 



9x9 


(3.57) 
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This was unappetizing to invert, but can be done using the methods developed in the 
previous sections. The inverse is 


$ = u — -9 x uj + gd x [6 x w) 

It 

= (1 — g\0\ 2 )u — x w + g{6 ■ u)$ 


(3.58) 


where 


_ sin(|^|/2) - (|0|/2)cos(|0|/2) 
9(|,|) |»Psin(|»|/2) 


(3.59) 


The matrix equivalent of (3.58) is 


°Z = (1 - <7| o 0| 2 )I 3 - ^[ a 0] + g a 9 a 0 J 


(3.60) 


where we have added the node designation a for clarity. 

Using the same arguments as before, we have for the internal forces the relation 


°f= a f+ - a 9 x a f+g a 0x ( a tf x a f). 

It 


(3.61) 


relating the old internal forces a f to the new “f. 

It is straightforward but very tedious to derive the second variation, as follows: 


Kmat = Z r KZ 

K seom = ~^[ a f]°Z + g{( a 0 • a f) a Z + a 0 ( a f a Z) - 2 a f ( a $ Ta Z)} 

+ * (** x af )> { a ° TaZ > 

where the derivative of g is expressed as 


(3.62) 


g ' + sin 0 cos ^) — 2 sin 2 ^ 

a 9\ 16^ 4 sin 2 4> 


(3.63) 


with the definition 

l a 0l 

<t>= -J- ( 3 * 64 ) 

Note that this linearization is more complicated than that involving the pseudovector with 
the tangent normalization, g ' requires special treatment for small angles, accomplished 
by using a power series expansion. All that was said before applies here, except for the 
additional complications in (3.61) through (3.64). 
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3.2.3. 2. Testing with an Example Potential 

We used the same potential as in the last section, but this time interpreting the pseudovec- 
tor as “angles” requiring the linearization just introduced. Here, in contrast to the last 
problem, there is a special case where the answer comes out in just one iteration. If one 
begins with the identity as the initial nodal triad (and some nonzero $o ), the increments 
produced solve this quadratic potential exactly. When these increments are turned into a 
surface triad, and that triad multiplies the starting value (the identity), we are left with 
something that, when inverted, reproduces the same angles that were solved for in the first 
place. Obviously, for this special case, one iteration is sufficient. On the other hand, any 
other starting position converts the linear problem in 0 to a nonlinear problem (because 
rotations neither add nor commute). This problem converges quadratically regardless of 
the starting state to the correct solution. Convergence was absolutely insensitive to any 
nonzero starting position. 

As before, the stiffness can be easily shown to be symmetric when equilibrium is 
satisfied, and, as before, the solution takes place at the same rate whether the stiffness is 
symmetrized or not. It is to be noted that the “super convergence” obtained by identity as 
the starting state requires a symmetrized stiffness, because only then does the contribution 
from Z vanish and the linear problem is retrieved. At this point, we decided to try a 
problem with fixed moments to see what effect this nonconservative external force had 
on the solution. In this case, of course, the skew-symmetric part of the stiffness does not 
vanish, but instead is proportional to one half the external moments. As expected, we 
obtained the rapid, quadratic convergence only when the unsymmetric stiffness was used. 
Otherwise, the convergence was linear at best. We note here in passing that unless 

m-60 x 6<l> = 0 (3.65) 

where m is the unbalanced moment and where 60, 6<f> are two independent allowable 
variations of the rotational freedoms, the forcing is nonconservative. This result is from 
[9], but it has appeared in the literature many times. From (3.65), any boundary condition 
that restrains only one rotation is not conserving, whereas a free edge and a boundary 
allowing only one free rotation are conservative. One can visualize why the boundary with 
the single rotational constraint is nonconserving with the following argument: 

1) Consider the freedom in the x direction as forbidden, i.e. no rotation about that 
axis can take place. 

2) Now rotate about the z axis 90° (allowed), followed by a rotation of some 0 about 
the y axis. Then rotate the result back by -90° about the z axis again. 

3) The result is a rotation by 0 about the forbidden x direction! 

Certainly, less extreme angles lead to similar results; this boundary is path-dependent and 
therefore nonconservative. ARGYRIS [1,2] has treated this subject at length, from which 
he derived his conserving semi-tangential rotations and symmetric geometric stiffnesses. 
The problem of nonconservative boundaries will come back to haunt us in succeeding 
sections. 
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3. 2.3.3. Implications for General Finite Element Analysis 

The arguments in subsection 3.2. 1.6 also apply here. The generality of what we have 
derived makes these results immediately applicable to the finite element code that uses 
elements that depend explicitly on pseudovectors with either the tangent normalization, 
or normalized as angles. As before, no modifications to the sequence of calculations are 
necessary, provided that 

1) The local element pseudovectors have normalization consistent with the choice of 
linearization from Section 3.2.2 or 3.2.3. 

2) The modifications to the linearization described above are implemented. Note that 
the band structure and coupling are exactly as in the case of the potential as a 
function of unit vectors. The additional computational effort is the same as before. 

3) The update formula for rotations is the exponential of a pseudovector normalized 
as angles, just as in the case of the unit vector potential. The modification we had 
for the C° elements works here, too. 

Whereas for the unit vector case (C° elements), the right hand side remains unchanged 
from what exists now, for elements that depend on angles like the 410, the right hand side 
is altered. This is, of course, unsettling. Fortunately, the grids required for reasonable 
convergence usually imply small local element angles. In that case, °Z approaches the 
identity, and the right hand side as modified differs little from the original. Convergence 
with respect to discretization remains intact. On the other hand, the corrections here cer- 
tainly enhance such convergence, allowing for coarser grids and hence potential substantial 
savings in computational and manpower costs. And, as we shall see, the faster, more stable 
convergence to the solution will yield further benefits, including the possibility of solving 
problems heretofore impossible because of poor convergence. 


3.3. Final Steps in the Linearization 

We are almost ready to put all this together into a working package, and demonstrate 
improved convergence for the corotational implementation. But first, since we know that 
the elements as now configured only approximately obey local equilibrium in the deformed 
configuration, we must come to terms with the effect of the projection operator P on the 
second variation. By utilizing the global invariance properties of the geometric stiffness, 
we will be able to account for all or almost all of the effect of the variation of P and the 
local element frame. 


3.3.1. The Behavior of the System under Rigid Rotation 

If the coordinate frame is rotated by a given infinitesimal angle -6u, internal and boundary 
forces of the system will appear to rotate by +6u. The relation, from Eq. (3.4) is 

6 a f =6uiX a f 
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But the stiffness matrix is the first order term in the Taylor expansion of the internal force, 
so that certainly 

K(<t> • w) = x “ x “f r ,-”} T 

= -{-,|f<r|,|T,|,-} r w (3.66) 

= — [F]w 

where the symbol [Fj denotes a long column vector consisting of a 3 x 3 skew-symmetric 
matrix for the translational a f t r and the rotational a f r internal forces for each node a in 
the system. Here the left superscripts denote node numbers, and is the infinitesimal 
rotation vector as defined in Section 2. We know that such a rotation is not conserving, 
because moments fixed in space also rotate, and because rotations do not commute. We 
shall see the implications of this in the next section. 

3.3.2. The End Run: the Completion of the Linearization 

Let us for the moment suppose we knew in advance the true tangent stiffness, and call this 

A 

one K. From the definition of the projector 

p = I - <h T 

we can “solve” for the identity 

I = P + <fa T (3.67) 

and pre- and postmultiply the stiffness to yield 

K = (P r + 70 T )K(P + <h T ) (3.68) 

Multiplying this out, we have 

K = P t KP + 70‘KP + P T K<h T + 1<t> T K<h T (3.69) 

We invoke (3.66) to yield 

K = P r KP - 7[F] t P - P t [F]7 T - 7[F \ T <h T (3.70) 

The first term on the right hand side can be identified as the “material” stiffness with the 
“true” stiffness replaced with the one we have calculated and projected with P as 

K = P r KP (3.71) 

where K results from ordinary finite element computations. The approximate version of 
(3.70) becomes 

K = K - 7(F1 t P - P r (Fh r - 7[F] T 07 T (3.72) 

Note that the added terms depend only on the internal forces. We started out with a 
symmetric matrix and assumed that by invoking the global rotation properties, we would 
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end up with a symmetric matrix. We see, however, that by reintroducing the definition of 
P we get instead 

K = K-7[F] r -[F]7 r + 70 T [F]7 r (3.73) 

which cannot be symmetric unless the 3x3 matrix 

Kred = * r [F] (3.74) 

is also symmetric. Using the properties (3.66) we find a contradiction: the matrix is 

unsymmetric, which means (3.66) describes an nonconservative system. We now take a 
cue from [1], and assume that the moments rotate half as much as the forces. This is what 
he calls the “semi-tangential” moment; his geometric stiffnesses rotate the forces by — |F], 
where 

= (3.75) 

It can be shown directly that K re d is symmetric when [F] replaces [F], provided the element 
internal forces are self-equilibrated. Since the projector P enforces local equilibrium with 
respect to the deformed configuration, this requirement is always satisfied, so that K re d is 
symmetric, as is the total tangent stiffness. 

To summarize the operations, we begin by calculating the ordinary, consistent tangent 
stiffness K, including the parts of the linearization derived in Sections 3.2 that apply to the 
particular situation (such as for C° elements as a function of a unit vector, for example). 
We also calculate the consistent internal force, also using the consistent linearization. Only 
then do we apply the projector to both the internal force and stiffness (see Eq. (3.71)). 
The internal forces are divided into two classes. For each set of translational freedoms, 
they are arranged into a skew-symmetric array, and placed into their proper position in 
[F]. Then half the rotational internal forces are likewise arranged and placed into their 
proper place in [Fj. We calculate the symmetric matrix 


K„i = * T [F] (3.76) 

and then the total tangent stiffness by 

K (arv = P T RP - -r|F) T - (F]i r + iK r , d i T (3.77) 


3.4. The Proof in the Pudding: Testing in the STAGS Code 

With all the pieces in place, we are in the position to implement the new theory in a full- 
blown finite element code. At this time, only one element, the often-used SH410 was tried. 
Software for the linearization of the element energy as a function of a pseudovector was 
extremely simple to implement, involving a very small amount of programming. Existing 
utilities carried us the bulk of the way. The operations to construct P and to implement 
(2.18) and the first term in (3.77) was in place from the linear work in Section 2, except 
that we had to pass deformed coordinates instead of undeformed nodal positions into the 
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various calculations for the projector. The remaining terms in (3.77) required only a little 
additional effort. 

Our first test was to produce the numerically-differentiated tangent stiffness and com- 
pare it with our calculated version. When we did, we got poor agreement! When we 
ran the same as a test case, the solution converged only linearly, even with true Newton 
selected as the solution option. The printed numerical stiffness was still unsymmetric, just 
as before. But right away we noticed that we had used as our external load consistent 
line loads derived from moderate rotation theory. These loads involve moments about 
fixed axes. What we have said about nonconservative systems certainly applies here: a 
quick check on the stiffness demonstrated that it was out of symmetry by exactly half the 
unsatisfied external moments. 

We repeated the runs with another problem with conserving boundaries and fixed 
point loads. Immediately, the numerically-calculated stiffness became symmetric. The 
calculated stiffness agreed to within 5 orders of magnitude (relative to the principal di- 
agonals) of the stiffness calculated by finite difference. Finally, and most important, we 
regained the true Newton quadratic convergence for the entire solution history, including 
when the structure was highly deformed and collapsed (past the limit point). Some small 
differences were seen between certain stiffness terms, a subject to be pursued later with 
the total derivative of the projection operator. We will then hope to show that the result 
will be almost, if not exactly, what we have derived here. One must note that unless all 
steps are included in the linearization, the good convergence is lost. This includes the 
modifications to the corotational software described in Sections 3. 2. 1.6 and 3. 2.3.3. For 
the STAGS 410 element, the linearization consistent with the corotational implementation 
described in [7] is with the pseudovector with the tangent normalization. This means, in 
particular, the use of relations (3.53), (3.54), (3.55), and taking note of the band structure 
(3.35), carrying out the operation in (3.28) with the value of a Z from (3.55). Any terms 
left out because of assumed “smallness” eventually destroy the quadratic convergence for 
solutions in the large deformation range. 

The next runs were made using the normalization of the pseudovector as angles. After 
linking to routines that produce the appropriate linearization, the convergence was just as 
good, and so were the answers. No real difference was observed between the two sets of 
runs. 

It is hard at this point to gauge what effect this will have on the ordinary, garden 
variety analysis. We have shown already that there is no profound difference in the over- 
all converged solution. This is good, because it means that the corotational theory as it 
stood was adequate for most situations. This fact has been buttressed by a multitude 
of runs that have been completed satisfactorily. On the other hand, having a consistent 
linearization is interesting in its own right. Even if benefits to the ordinary user are not 
dramatic, confidence that the solution is consistently derived is comforting, especially since 
the right hand side is also modified, both through the projector P and the linearization 
of the pseudovector (for elements like 410). One can certainly foresee cases where conver- 
gence previously was difficult, and therefore solution impossible with things the old way. 
In particular, when taking advantage of the new Thurston Equivalance Transformation 
processor (ET) [11] to carry a solution through bifurcation in the presence of modal inter- 
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action, a consistent tangent stiffness can be crucial. The modifications to the right hand 
side from the projector or for functions of a pseudovector required for the SH410 element 
will also be important when calculating the higher-order terms needed to properly account 
for the modal interaction; the combined effect of an inconsistent tangent and approximate 
higher-order terms could put the solution in jeopardy. We anticipate experimenting with 
collapse and bifurcation problems that were difficult in the past to see if consistency will 
gain us the robustness to overcome the complicated deformations characteristic of these 
structures. 

To this point, we have demonstrated that even for the modified Newton solution 
strategy, we get a ratio of about .85 in the iteration count as compared to 1 for the old 
way when just the projector is introduced. Since element-related computation accounts 
for the bulk of the work for these simple problems, there is a net loss in overall efficiency, 
with a ratio of about 1.1. When the very inexpensive additions to the linearization are 
added, the ratio of iterations drops to about .65, with the work level to about .95 for some 
simple cases. If the iteration ratio holds up, for large problems where matrix operations 
tend to dominate (none of these were altered in any way), we may approach a gain of .65 
to one in efficiency, measured as relative time spent doing the calculations. These results 
were computed using a simple, but nontrivial, small problem. For those problems that 
are intrinsically difficult, the difference may be between success and failure. In the case of 
true Newton, its status has been restored as a robust method of solution for conservative 
problems defined by an elastic potential. Before, true Newton cost may times more than 
modified Newton. For this small case, the factor was reduced to about 1.5, which of course 
will be much greater for large problems; true Newton is always expensive because we must 
factor every iteration. 

3.5. Extending the Conservative Boundary Condition 

The one messy aspect of things as they now stand is the possibility of the user generating 
a nonconservative forcing or boundary condition without knowing it. Certainly, if one 
has true follower forces or localized nonconserving pressure loads, one cannot expect to 
avoid the consequences. For pressure loads derived from a potential, we are in good shape, 
because a load stiffness is introduced to include the effects of the change in direction of 
the pressure load. The implementation in STAGS is not affected by what we have derived 
here. 

For other types of load, we have a problem. We have already discussed how trouble- 
some boundaries can be. Loads determined by moments fixed in space caused delays in 
the research presented here. Fortunately, the solution to these problems is already implied 
by Eq. (3.56). There is a term in that equation of the form 

$ T Ke 0 = 0 T m (3.78) 

where certainly m could be interpreted as an external moment. But (3.78) is a potential, a 
function of a well-defined, unique pseudovector with, in this case, the angle normalization. 
Therefore, (3.61) through (3.64) apply. What one gets is a modification to the applied 
load, as expected. In addition, however, we get the new load stiffness terms (3.62). All 
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angles of course would be based on the relative rotation the affected node has undergone 
from its resting state ( a S(i)). One only has to show that the external load defined in 
this way is consistent with those cases that satisfy (3.65). But (3.65) states basically the 
moment must follow the direction of the rotational freedom variation. In practice, this 
means that either the moment is zero, or that rotation is possible only about a single axis. 
If the moment is zero, (3.62) contributes nothing. If the rotation is permitted to take place 
about a single axis, one can immediately see by examining (3.58) that 


0 x u> = 0 


because only one direction is permitted, and hence 


= uj 


(3.79) 


Therefore, no further transformations are needed, and the usual interpretation of the ex- 
ternal load survives. Note that because of the scalar terms in (3.53), such correspondence 
is impossible when the tangent pseudovector is selected. 


For follower nodal loads, we have available only the contribution to the residual, as 
follows: 



a f r = 0 


(3.80) 


where °f is a given nodal follower load rigidly attached to the nodal triad. The load stiffness 
is derived directly from (3.79), (3.4), and (2.13): 


K 


gtom 


o — [ a f] 

0 0 


(3.81) 


This matrix is unsymmetric for all cases because the follower load is nonconservative. 


3.5.1. Application to Arc-Length Solution Algorithms 

One important constraint condition to be considered arises whenever continuation param- 
eter solution procedures are invoked. In these procedures, the load parameter is cast as 
dependent variable to be solved for along with the displacements. Typically, the additional 
constraint equation is a function of some “arc length” along an approximate tangent to 
the solution path, often approximated by the relation 

t r (u (i) - u (0) ) + *a(A (0 - A (0) ) - At; = 0 (3.82) 

where U(,) are the displacements of the current iterate *, where (0) refers to the latest 
converged solution ( iteration 0), and where t is an approximate “tangent” to the current 
solution path. Here, A is the load factor and t\ is that component of t belonging to the 
load factor, and A 77 is the arc step length (independent variable) selected beforehand. 
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The problem with (3.82) is that the difference between the displacement states is not 
valid for rotational freedoms when the arc step is large, as can be the case for a fine-tuned 
algorithm. Fortunately, we can bring the results of the previous sections to bear on the 
constraint equation if we redefine (3.82) as follows: 


t t^( u (») ~ u (o)) + f^(R) + *a (A( i) - A( 0 )) - At/ = 0 (3.83) 

where 

R = ®(‘)®(o) (3.84) 

is the relative rotation from the reference state 0 to the current iterate i. Note that again 
we have redefined the system state subscripts to refer to iteration count , with 0 as the 
reference converged solution. If we define 0 by the relation 

exp([0]) = R (3.85) 

then (3.83) reverts to (3.82) for small angles. Following the arguments in Sect. 3.2.2, we 
see that Eq. (3.44) or (3.57) applies, where this time 6 is the appropriate pseudovector 
corresponding to the relative rotation R through Eq. (3.85). 

One must evaluate (3.83) and its first derivative to complete the arc-length calcula- 
tions. To evaluate (3.83), the relative rotation is calculated from (3.84) and used to obtain 
0 (3.85) by standard means. Either the angle or tangent normalization can be used. The 
derivative of (3.83) is treated just like the right hand side of (3.78), from which a modified 
“internal force” is extracted through the use of (3.53) as follows: 

“t m = a t + x a t + \{ a 0- a t) a 0 (3.86) 

for the pseudovector 6 with tangent normalization corresponding to 0 in (3.85). For angle 
normalization, we proceed using (3.61) instead. In either case, “t m replaces t in the arc- 
length extended stiffness for the rotational freedoms pertaining to node a. The rest of the 
calculations are unchanged. 
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4. Summary and Conclusions 

In this work, we have succeeded in largely completing the element-independent corotational 
formulation introduced in [7]. We have demonstrated that the theory as it has existed up 
to now is incomplete. We have uncovered these deficiencies and derived what is necessary 
to fill the gaps in the theory. Sample problems were introduced that cover the relevant 
points, and what was derived was tested on these problems. The linearization of a potential 
as a function of a unit vector, or of a pseudovector with two kinds of normalization was 
shown to be completely consistent, with true Newton, quadratic convergence to solution. 
Finally, the parts were put together in the large structural analysis code STAGS and 
superior convergence was demonstrated. Extensions to follower loads and unconservative 
boundaries were derived. 

This work culminates a long effort to realize the full potential of the element- 
independent corotational theory. As originally posed, such a theory would greatly pro- 
mote a host of inexpensive finite elements to robust performance for difficult nonlinear 
collapse problems. The roadblocks have been fundamental element inconsistencies (such 
as nonsastisfaction of equilibrium) and the accompanying poor convergence properties. 
Two important gaps have now been filled: 

1) The element is automatically self- equilibrated, even if in its original conception local 
equilibrium was not satisfied. 

2) The linearization is consistent. 

What we have not done is fill in such holes as rank-deficiency, distortion sensitivity in a 
plane, or shear locking. But we envision that perhaps elements could be defined without 
concern for frame invariance, something to be enforced by the methods developed here. 
Perhaps trigonometric functions could once again be used as shape functions for shell 
elements? Finally, if it could be possible to define the element interpolation entirely as a 
function of scalars defined by inner products of the unit nodal normals and the displacement 
vectors taken in proper combination, frame invariance would be guaranteed, and no added 
“corotational” theory would be necessary. What we have here provides the recipe for the 
linearization of these elements also. 
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5. Suggestions for Future Research 

The most immediate need is to test thoroughly what we have already in place. We wish 
to see if more difficult problems can now be tackled. Following this, careful benchmarks 
should be run to measure what immediate gains in solution efficiency have been realized. 
In addition, the sensitivity of solutions to the order, sequence, and frequency of the element 
frame update should be investigated to see whether moderate rotation theory within the 
local element environment is in effect. If it proves true that the solution is insensitive to 
the frequency of the element frame update, strategies for taking advantage of this should 
be investigated. Provision for the two normal rotations in the special SH411 elements 
should be made before any of this work is implemented for that element. 

Before publishing the complete linearization including the derivative of projection 
operator, we wish to calculate all derivatives analytically, using symbolic manipulation on 
the computer. Following the completion of that effort, attention should be paid to the 
nonconserving boundaries and external loads. 

Areas of new research abound. For example, we mentioned elements as a function 
of fundamental relative scalar quantities that are already frame-invariant. Trigonometric 
shape functions could be exhumed, and their use as shape functions should be examined 
in the light of what we can do now. Elements completely lacking invariance to rotation or 
translation could be used with the generalized projection operator (the translational part 
is trivial) to exploit the otherwise favorable behavior they might have. 
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Normalized maximum displacement 



Pinched Hemisphere Problem 

R/h ■ 250 



Figure 1. 


Pinched hemisphere results for linear and nonlinear analysis with SH410 
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Normalized maximum displacement 


Pinched Hemisphere Problem 

R/h ■ 250 



Figure 2. Pinched hemisphere results for various elements with and without P 
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Figure 3. Pinched cylinder results for various elements with and without P 
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